function [theta_est1, theta_est2, se, p_value]=Master_Fox_mixed_discrete_elecavg_3param_grID_wtp_j_bs(subsample,inc,set_seed,matlab_seed,get_nb_draws, j, index_bs)
																										 
																									
tic

global state week J T S MT MTvector K WB param_nopt...
price_num rebate_estar_num kwh_num elec_num estar_num standard0_num standard10_num standard15_num standard20_num standard25_num rank_num ...
price_denom rebate_estar_denom rebate_denom kwh_denom elec_denom estar_denom standard0_denom standard10_denom standard15_denom standard20_denom standard25_denom rank_denom ...
bchoice bchoice_s bchoice_ij id_i id_j id_ij id_m id_n  demo pidnest choicenest_ij brandweek_denom brandweek_num prod_denom prod_num denom_tmp2 denom_tmp3 denom...
type_id2_3_demo1_num type_id2_3_demo1_denom type_id3_demo1_num type_id3_demo1_denom AV_demo1_num AV_demo1_denom  ice_sc_demo1_num ice_sc_demo1_denom...
type_id2_3_demo2_num type_id2_3_demo2_denom type_id3_demo2_num type_id3_demo2_denom AV_demo2_num AV_demo2_denom  ice_sc_demo2_num ice_sc_demo2_denom...
type_id2_3_demo3_num type_id2_3_demo3_denom type_id3_demo3_num type_id3_demo3_denom AV_demo3_num AV_demo3_denom  ice_sc_demo3_num ice_sc_demo3_denom...
type_id2_3_demo4_num type_id2_3_demo4_denom type_id3_demo4_num type_id3_demo4_denom AV_demo4_num AV_demo4_denom  ice_sc_demo4_num ice_sc_demo4_denom...
type_id2_3_demo5_num type_id2_3_demo5_denom type_id3_demo5_num type_id3_demo5_denom AV_demo5_num AV_demo5_denom  ice_sc_demo5_num ice_sc_demo5_denom...
type_id2_3_demo6_num type_id2_3_demo6_denom type_id3_demo6_num type_id3_demo6_denom AV_demo6_num AV_demo6_denom  ice_sc_demo6_num ice_sc_demo6_denom;
global id;

global LLC_5_18;

nb_draws=get_nb_draws;

date_now = clock;
date_now = strcat(num2str(date_now(1)),'_',num2str(date_now(2)),'_', num2str(date_now(3)), num2str(date_now(4)), num2str(date_now(5)));
%diary(strcat('\\c3\rdat\SHoude\Research/sears/estar_results/struct_est/','log_','Fox_discrete_FEdemoshrt_WTP_tau1_',date_now,'.log'));

%Create Output: See below and make changes
results_file1=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Results_Estimation\','matrix_Fox_discrete_bup_mixedFEdemoshrt_WTP_v',num2str(subsample),'_inc',num2str(inc),'_sd',num2str(set_seed),'_cty_tau1_grid_v3_','bs_',num2str(j),'.mat');
results_file2=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Results_Estimation\','Fox_discrete_bup_mixedFEdemoshrt_WTP_v',num2str(subsample),'_inc',num2str(inc),'_sd',num2str(set_seed),'_cty_tau1_grid_v3_','bs_',num2str(j),'.mat');

%Read Files
	tmp_file1=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\choiceset_identifier_trimester_week_zipcode_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    id=load(tmp_file1);
    id=id(index_bs,:);
    id(:,1:7)=[];
    choicest_size=sum(id,2);
    %id_discard=find(choicest_size>200 | choicest_size<50);
    id_discard=find(choicest_size>400);
    id(id_discard,:)=[];  
    choicest_size(id_discard)=[]; 
    id=id';
    [id_i,id_j,id_s]=find(id);
    id_ij=find(id);
    [id_m,id_n] = size(id);
    
    
 	%bchoice=load('H:\Research\sears\estar_data\refrigerators\bchoice_046_2008_2012_v11022017_struct_52171_4.csv'); 
    tmp_file2=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\bchoice_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
 	bchoice=load(tmp_file2);
 	bchoice=bchoice(index_bs,:); 
    bchoice(:,1)=[];
    bchoice(id_discard,:)=[]; 
    bchoice=bchoice';
    [bchoice_i,bchoice_j,bchoice_s]=find(bchoice);
    bchoice_ij=find(bchoice);
    [bchoice_m,bchoice_n] = size(bchoice);  
    
    bchoiceIJ=bchoice(id_ij);
    MT=size(id,2);
    MTvector=ones(MT,1);
    J=size(id,1);
    
    tmp_file3=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\estar_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    estar=load(tmp_file3); 
    estar=estar(index_bs,:);
    estar(:,1:7)=[];
    estar(id_discard,:)=[]; 
    estar=estar';
    estar_num=estar(bchoice_ij);
    estar_denom=estar(id_ij);
    clear estar;
    
    tmp_file4=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\rebate_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    rebate=load(tmp_file4); 
    rebate=rebate(index_bs,:); 
    rebate(:,1:7)=[];
    rebate=rebate/1000;
    rebate_J=repmat(rebate,1,J);
    rebate_J(id_discard,:)=[]; 
    rebate_J=rebate_J';
    rebate_denom=rebate_J(id_ij);
    rebate_estar_num=rebate_J(bchoice_ij).*estar_num;
    rebate_estar_denom=rebate_J(id_ij).*estar_denom;
    clear rebate rebate_J;
    
    tmp_file5=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\tax_rate_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tax=load(tmp_file5); 
    tax=tax(index_bs,:);
    tax(:,1:7)=[];
    tax(id_discard,:)=[]; 
    taxp=1+repmat(tax',J,1);
    
    tmp_file6=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\price_mean_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    price_net=load(tmp_file6); 
    price_net=price_net(index_bs,:);
    price_net(:,1:7)=[];
    price_net(id_discard,:)=[]; 
    price_net=price_net';
    price=taxp.*price_net;
    price=price/1000;
    price_num=price(bchoice_ij);
    price_denom=price(id_ij);
    clear price_net taxp tax;
    
    tmp_file7=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\kwh_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    kwh=load(tmp_file7); 
    kwh=kwh(index_bs,:);
    kwh(:,1:7)=[];
    kwh(id_discard,:)=[]; 
    kwh=kwh';
    kwh=kwh/1000;
    kwh_num=kwh(bchoice_ij);
    kwh_denom=kwh(id_ij);
    clear kwh;
    
    tmp_file8=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\elec_county_046_avg_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    elec=load(tmp_file8); 
    elec=elec(index_bs,:);
    elec(:,1:7)=[];
    elec_J=repmat(elec,1,J);
    elec_J(id_discard,:)=[]; 
    elec_J=elec_J';
    elec_num=elec_J(bchoice_ij).*kwh_num;
    elec_denom=elec_J(id_ij).*kwh_denom;
    clear elec_J;
    
    tmp_file9=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\demo_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    demo=load(tmp_file9); 
    demo=demo(index_bs,:);

%hd_id income2 income3 education2 education3 fam_size age political2 political3
    demo(:,1)=[];
    demo(:,1:2)=[];
    demo(:,3)=demo(:,3)/10;
    demo(:,4)=demo(:,4)/100;
    demo1_J=repmat(demo(:,1),1,J); %education2
    demo2_J=repmat(demo(:,2),1,J); %education3
    demo3_J=repmat(demo(:,3),1,J); %fam_size
    demo4_J=repmat(demo(:,4),1,J); %age
    demo5_J=repmat(demo(:,5),1,J); %political2
    demo6_J=repmat(demo(:,6),1,J); %political3
   
%Create Interactions Attributes-Demographics     
	tmp_file10=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\type_id_2_3_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
	type_id2_3=load(tmp_file10); 
	type_id2_3=type_id2_3(index_bs,:);
    type_id2_3(:,1:7)=[];
   
    type_id2_3_demo1=type_id2_3.*demo1_J;
    type_id2_3_demo1(id_discard,:)=[]; 
    type_id2_3_demo1=type_id2_3_demo1';
    type_id2_3_demo1_num=type_id2_3_demo1(bchoice_ij);
    type_id2_3_demo1_denom=type_id2_3_demo1(id_ij);
    clear type_id2_3_demo1;
    
    type_id2_3_demo2=type_id2_3.*demo2_J;
    type_id2_3_demo2(id_discard,:)=[];
    type_id2_3_demo2=type_id2_3_demo2';
    type_id2_3_demo2_num=type_id2_3_demo2(bchoice_ij);
    type_id2_3_demo2_denom=type_id2_3_demo2(id_ij);
    clear type_id2_3_demo2;
    
    type_id2_3_demo3=type_id2_3.*demo3_J;
    type_id2_3_demo3(id_discard,:)=[];
    type_id2_3_demo3=type_id2_3_demo3';
    type_id2_3_demo3_num=type_id2_3_demo3(bchoice_ij);
    type_id2_3_demo3_denom=type_id2_3_demo3(id_ij);
    clear type_id2_3_demo3;
    
    type_id2_3_demo4=type_id2_3.*demo4_J;
    type_id2_3_demo4(id_discard,:)=[];
    type_id2_3_demo4=type_id2_3_demo4';
	type_id2_3_demo4_num=type_id2_3_demo4(bchoice_ij);
    type_id2_3_demo4_denom=type_id2_3_demo4(id_ij);
    clear type_id2_3_demo4;
    
    type_id2_3_demo5=type_id2_3.*demo5_J;
    type_id2_3_demo5(id_discard,:)=[];
    type_id2_3_demo5=type_id2_3_demo5';
    type_id2_3_demo5_num=type_id2_3_demo5(bchoice_ij);
    type_id2_3_demo5_denom=type_id2_3_demo5(id_ij);
    clear type_id2_3_demo5;
    
    type_id2_3_demo6=type_id2_3.*demo6_J;
    type_id2_3_demo6(id_discard,:)=[];
    type_id2_3_demo6=type_id2_3_demo6';
    type_id2_3_demo6_num=type_id2_3_demo6(bchoice_ij);
    type_id2_3_demo6_denom=type_id2_3_demo6(id_ij);
    clear type_id2_3_demo6;
   
   
    tmp_file11=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\AV_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
    AV=load(tmp_file11);
    AV=AV(index_bs,:); 
    AV(:,1:7)=[];
    AV=AV/100;
    
    AV_demo1=AV.*demo1_J;
    AV_demo1(id_discard,:)=[];
    AV_demo1=AV_demo1';
    AV_demo1_num=AV_demo1(bchoice_ij);
    AV_demo1_denom=AV_demo1(id_ij);
    clear AV_demo1;
    
    AV_demo2=AV.*demo2_J;
    AV_demo2(id_discard,:)=[];
    AV_demo2=AV_demo2';
    AV_demo2_R=repmat(AV_demo2,1,nb_draws);
    AV_demo2_num=AV_demo2(bchoice_ij);
    AV_demo2_denom=AV_demo2(id_ij);
    clear AV_demo2;
    
    AV_demo3=AV.*demo3_J;
    AV_demo3(id_discard,:)=[];
    AV_demo3=AV_demo3';
    AV_demo3_num=AV_demo3(bchoice_ij);
    AV_demo3_denom=AV_demo3(id_ij);
    clear AV_demo3;
    
    AV_demo4=AV.*demo4_J;
    AV_demo4(id_discard,:)=[];
    AV_demo4=AV_demo4';
    AV_demo4_num=AV_demo4(bchoice_ij);
    AV_demo4_denom=AV_demo4(id_ij);
    clear AV_demo4;
    
    AV_demo5=AV.*demo5_J;
    AV_demo5(id_discard,:)=[];
    AV_demo5=AV_demo5';
    AV_demo5_num=AV_demo5(bchoice_ij);
    AV_demo5_denom=AV_demo5(id_ij);
    clear AV_demo5;
    
    AV_demo6=AV.*demo6_J;
    AV_demo6(id_discard,:)=[];
    AV_demo6=AV_demo6';
    AV_demo6_num=AV_demo6(bchoice_ij);
    AV_demo6_denom=AV_demo6(id_ij);
    clear AV_demo6;
   
    clear AV;
    
    tmp_file12=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\ice_sc_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
    ice_sc=load(tmp_file12); 
    ice_sc=ice_sc(index_bs,:); 
    ice_sc(:,1:7)=[];
   
    ice_sc_demo1=ice_sc.*demo1_J;
    ice_sc_demo1(id_discard,:)=[];
    ice_sc_demo1=ice_sc_demo1';
    ice_sc_demo1_num=ice_sc_demo1(bchoice_ij);
    ice_sc_demo1_denom=ice_sc_demo1(id_ij);
    clear ice_sc_demo1;
    
    ice_sc_demo2=ice_sc.*demo2_J;
    ice_sc_demo2(id_discard,:)=[];
    ice_sc_demo2=ice_sc_demo2';
    ice_sc_demo2_num=ice_sc_demo2(bchoice_ij);
    ice_sc_demo2_denom=ice_sc_demo2(id_ij);
    clear ice_sc_demo2;
    
    ice_sc_demo3=ice_sc.*demo3_J;
    ice_sc_demo3(id_discard,:)=[];
    ice_sc_demo3=ice_sc_demo3';
    ice_sc_demo3_num=ice_sc_demo3(bchoice_ij);
    ice_sc_demo3_denom=ice_sc_demo3(id_ij);
    clear ice_sc_demo3;
    
    ice_sc_demo4=ice_sc.*demo4_J;
    ice_sc_demo4(id_discard,:)=[];
    ice_sc_demo4=ice_sc_demo4';
    ice_sc_demo4_num=ice_sc_demo4(bchoice_ij);
    ice_sc_demo4_denom=ice_sc_demo4(id_ij);
    clear ice_sc_demo4;
    
    ice_sc_demo5=ice_sc.*demo5_J;
    ice_sc_demo5(id_discard,:)=[];
    ice_sc_demo5=ice_sc_demo5';
    ice_sc_demo5_num=ice_sc_demo5(bchoice_ij);
    ice_sc_demo5_denom=ice_sc_demo5(id_ij);
    clear ice_sc_demo5;
    
    ice_sc_demo6=ice_sc.*demo6_J;
    ice_sc_demo6(id_discard,:)=[];
    ice_sc_demo6=ice_sc_demo6';
    ice_sc_demo6_num=ice_sc_demo6(bchoice_ij);
    ice_sc_demo6_denom=ice_sc_demo6(id_ij);
    clear ice_sc_demo6;
   
    clear ice_sc;
    
      
    prod=repmat([1:J],MT,1)';
    prod_denom=prod(id_ij);
    prod_num=prod(bchoice_ij);
    
 clear data;   
toc  
%% 

 rho_5=1/(1+0.05);
 LLC_5_18=rho_5*(1-rho_5^18)/(1-rho_5);

%result_homFE=load('\\c3\rdat\SHoude\Research\sears\estar_results\struct_est\hom_FEdemoshrt_66000_allinc_hdonly_101.mat');
%result_homFE0=result_homFE.theta_est1;
%param_nopt=result_homFE0(1:J-1,1);
%hom_FEdemoshrt_wtp_LLC_5_18(11000)_inc(3)_sd(5)_05162017
% result_homFE=strcat('\\c3\rdat\SHoude\Research\sears\estar_results\struct_est\hom_FEdemoshrt_wtp_LLC_5_18',num2str(subsample),'_inc',num2str(inc),'_sd',num2str(set_seed),'_05162017.mat');
% eval(['load ' result_homFE]);
% result_homFE0=theta_est1;
% param_nopt=result_homFE0(1:J-1,1);
% 
% eta=result_homFE0(J)  
% tau=result_homFE0(J+1)
% pi=result_homFE0(J+2)
% theta=result_homFE0(J+3)
% beta=result_homFE0(J+4:J+15)

result_mixed=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Results_Estimation\','mixed_logit_FEdemoshrt_grID_wtp_100iters_',num2str(subsample),'_',num2str(inc),'_',num2str(set_seed),'_matlab_seed_',num2str(matlab_seed),'_nb_draws_',num2str(nb_draws),'bs_',num2str(j),'.mat');
eval(['load ' result_mixed]);
result_mixed0=theta_est1;
param_nopt=result_mixed0(1:J-1,1);

eta=result_mixed0(J)  
tau=result_mixed0(J+1)
pi=result_mixed0(J+2)
theta=result_mixed0(J+3)
beta=result_mixed0(J+4:J+15)


%%
  param_a(1,1)=eta;
  param_a(2,1)=tau;
  param_a(3,1)=pi;
  param_a(4,1)=theta;
  param_a(4+1:4+12,1)=beta;
  
  %grid_a=[-3:21]/7;
  %grid_a=[-2:6]/2;
  %eta_v=grid_a*eta-0.1;
  %theta_v=grid_a*theta-0.05;
  %theta_v=[-100,-80,-50,-33,theta,-25,-18,-12,-6,-1,-0.1,10,18];
  
  %Grid v1
  %theta_v=[-2,-1.5,-1.25,-1,-0.75,-0.5,-0.25,-0.1,0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.25,1.5,2,2.5,3,3.5,4];
  %eta_v=[-25,-18,-16,-14,-10,-7,-5,-4.5,-4,-3.5,-3.25,-3,-2.5,-2,-1,-0.1,1,2];
  
  %Grid v2	
  %theta_v=[-6,-5,-4,-3,-2.5,-2,-1.5,-1.25,-1,-0.75,-0.5,-0.25,-0.1,0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.25,1.5,2,2.5,3,3.5,4,5,6,7,8,9,10];
  %eta_v=[-25,-18,-16,-14,-10,-7,-5,-4.5,-4,-3.5,-3.25,-3,-2.5,-2,-1,-0.1,1,2];
  
  %Grid v3
  theta_v=[-3.5,-2.5,-2,-1.5,-1.25,-1,-0.75,-0.5,-0.25,-0.1,0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.25,1.5,2,2.5,3,3.5,4,5];
  eta_v=[-35,-30,-25,-18,-16,-14,-10,-7,-5,-4.5,-4,-3.5,-3.25,-3,-2.5,-2,-1,-0.1,1,2,3];
  
  
  %tau_v=[8*tau 6*tau 5*tau 4*tau 2.5*tau 1.75*tau 1.25*tau tau 0.75*tau 0.5*tau 0.25*tau 0 -0.25*tau -0.5*tau];
  %tau_v=[0*tau tau 2.5*tau];
  tau_v=tau;
  %grid_a=[-20:60]/10;
  %eta_v=grid_a*eta;
  %theta_v=grid_a*theta;
  %tau_v=[8*tau 6*tau 5*tau 4*tau 2.5*tau 1.75*tau 1.25*tau tau 0.75*tau 0.5*tau 0.25*tau 0 -0.25*tau -0.5*tau];
  
c=1
  for k=1:length(eta_v)
   for m=1:length(theta_v)
    for n=1:length(tau_v)
        
        grid_param(c,:)=[eta_v(k),theta_v(m),tau_v(n)];
   		c=c+1;     
    end
   end
  end

% ntasks=99;
dim_all=length(eta_v)*length(theta_v)*length(tau_v);
dim_all
for ii=1:dim_all
l=ii;   
	         eta_s=grid_param(l,1);
	         tau_s=grid_param(l,3);
	         theta_s=grid_param(l,2);
	         
	         param_a(1,1)=eta_s;    
	         param_a(2,1)=tau_s;    
	         param_a(4,1)=theta_s;
	            
	           ChoiceProbIJR(:,l)=ChoiceProbIJ_hom_FEdemoshrt_WTP_vs(param_a);    
                   
end

for ii=1:dim_all
l=ii   
         
         r18=SolveGeoSeries_shell(grid_param(l,1),grid_param(l,2));
         
         param_v(l,:)=[grid_param(l,1),grid_param(l,2),grid_param(l,3),r18];
       
         
end


eval(['save ' results_file1 ' ChoiceProbIJR param_v bchoiceIJ']); 

%%
%Estimation
tic

options_lsqlin=optimset('Display','iter','TolFun',1e-7,'MaxFunEvals',100000,'MaxIter',100000);
beta = lsqlin(ChoiceProbIJR,bchoiceIJ,'','',ones(1,dim_all),1,zeros(dim_all,1),ones(1,dim_all),rand(1,dim_all),options_lsqlin);
toc

%%
   
eval(['save ' results_file2 ' param_v beta']);    


end 






      





